module era5_regrid_mod

  use const_mod
  use latlon_data_mod
  use latlon_interp_mod
  use vert_interp_mod
  use era5_reader_mod

  implicit none

  private

  public era5_regrid_run

contains

  subroutine era5_regrid_run(model_data)

    type(latlon_data_type), intent(inout) :: model_data

    real(r8), allocatable, dimension(:,:  ) :: p0, t0
    real(r8), allocatable, dimension(:,:,:) :: u1, v1, t1, q1
    real(r8) t0_p, z0
    real(r8), parameter :: lapse = 0.006_r8 ! K m-1
    real(r8), parameter :: kappa = Rd / g
    real(r8), parameter :: lapse_kappa = lapse * kappa
    integer i, j, k, is, ie, js, je

    is = proc%local_mesh%lon_lb
    ie = proc%local_mesh%lon_ub
    js = proc%local_mesh%lat_lb
    je = proc%local_mesh%lat_ub

    allocate(p0(is:ie,js:je))
    allocate(t0(is:ie,js:je))
    allocate(u1(is:ie,js:je,num_era5_lev))
    allocate(v1(is:ie,js:je,num_era5_lev))
    allocate(t1(is:ie,js:je,num_era5_lev))
    allocate(q1(is:ie,js:je,num_era5_lev))

    call latlon_interp_bilinear(era5_lon, era5_lat, era5_mslp, proc%local_mesh, p0)
    call latlon_interp_bilinear(era5_lon, era5_lat, era5_t(:,:,num_era5_lev), proc%local_mesh, t0)
    t0_p = era5_lev(num_era5_lev)
    z0 = 0.0_r8
    do j = proc%local_mesh%lat_ibeg, proc%local_mesh%lat_iend
      do i = proc%local_mesh%lon_ibeg, proc%local_mesh%lon_iend
        t0(i,j) = t0(i,j) * (p0(i,j) / t0_p)**lapse_kappa
        model_data%ps(i,j) = p0(i,j) * (1.0_r8 - lapse * (model_data%zs(i,j) - z0) / t0(i,j))**(1.0_r8 / lapse_kappa)
      end do
    end do

    do k = 1, num_era5_lev
      call latlon_interp_bilinear(era5_lon, era5_lat, era5_u(:,:,k), proc%local_mesh, u1(:,:,k))
      call latlon_interp_bilinear(era5_lon, era5_lat, era5_v(:,:,k), proc%local_mesh, v1(:,:,k))
      call latlon_interp_bilinear(era5_lon, era5_lat, era5_t(:,:,k), proc%local_mesh, t1(:,:,k))
      call latlon_interp_bilinear(era5_lon, era5_lat, era5_q(:,:,k), proc%local_mesh, q1(:,:,k))
    end do

    do j = proc%local_mesh%lat_ibeg, proc%local_mesh%lat_iend
      do i = proc%local_mesh%lon_ibeg, proc%local_mesh%lon_iend
        call vert_interp_linear(era5_lev, u1(i,j,:), proc%local_mesh%lev, model_data%u(i,j,:), allow_extrap=.true.)
        call vert_interp_linear(era5_lev, v1(i,j,:), proc%local_mesh%lev, model_data%v(i,j,:), allow_extrap=.true.)
        call vert_interp_linear(era5_lev, q1(i,j,:), proc%local_mesh%lev, model_data%q(i,j,:), allow_extrap=.true.)
        call vert_interp_log_linear(era5_lev, t1(i,j,:), proc%local_mesh%lev, model_data%t(i,j,:), allow_extrap=.true.)
      end do
    end do

    deallocate(p0)
    deallocate(t0)
    deallocate(u1)
    deallocate(v1)
    deallocate(t1)
    deallocate(q1)

  end subroutine era5_regrid_run

end module era5_regrid_mod
